%--------------------------------------------------------------------------
% SIMULATION
%--------------------------------------------------------------------------

rng('default')
Nsim = 100; % number of firms


% simulate initial firm distribution without disasters,
% mimicking the 1947--2003 pre-estimation sample,
% based on which to initialize the actual simulation
burn_in         = 1000; % number of initial periods to discart
Tsim            = 10000;
[~,states_init] = sim_markov_chain_fast(P_nodis,Tsim+burn_in,1,1,1:2);
epsC_init       = randn(Tsim+burn_in,1);
epsCE_init      = [epsC_init,randn(size(epsC_init))]*chol([1,rho_CE;rho_CE,1]);
epsE_init       = epsCE_init(:,2);
g_sim           = mu_X(states_init) + sig_XS(states_init).*epsE_init + sig_Xid(states_init).*randn(Tsim+burn_in,Nsim); % [Tsim,Nsim]
k_sim_init      = NaN(size(g_sim));
k_sim_init(1,:) = k_bar(states_init(1));
k_QML_init      = NaN(size(g_sim));
k_QML_init(1,:) = k_bar(states_init(1));
last_issuance_init = ones(size(g_sim)); % times initial state of the Markov chain
for t = 1 : Tsim+burn_in-1 % maturity
    
    k_sim_init(t+1,:) = k_sim_init(t,:).*exp(g_sim(t,:));
    k_QML_init(t+1,:) = k_QML_init(t,:);
    
    default = k_sim_init(t,:) <= k_def(states_init(t));
    k_sim_init(t+1,default) = k_bar(states_init(t)).*exp(g_sim(t,default));
    k_QML_init(t+1,default) = k_bar(states_init(t));
    
    issuance = k_sim_init(t,:) >= k_iss(states_init(t));
    k_sim_init(t+1,issuance) = k_bar(states_init(t)).*exp(g_sim(t,issuance));
    k_QML_init(t+1,issuance) = k_bar(states_init(t));
    
    last_issuance_init(t:end,issuance)  = ones(Tsim+burn_in-t+1,sum(issuance)).*states_init(t);
    last_issuance_init(t:end,default)  = ones(Tsim+burn_in-t+1,sum(default)).*states_init(t);
    
end
k_sim_init(1:burn_in,:) = [];
k_QML_init(1:burn_in,:) = [];
last_issuance_init(1:burn_in,:) = [];
states_init(1:burn_in) = [];

% randomize cross-sections
ix = randperm(Tsim)';
k_sim_init = k_sim_init(ix,:);
k_QML_init = k_QML_init(ix,:);
last_issuance_init = last_issuance_init(ix,:);
states_init = states_init(ix);

% find periods corresponding to states 1 & 2
i1_init = find(states_init==1);
i2_init = find(states_init==2);


% simulate estimation sample
load aggregate_path_samples
delC = delC(:,1:5000);
epsC = epsC(:,1:5000);
states = states(:,1:5000);

Tsim   = size(states,1); % number of periods
Ssim   = size(states,2); % number of panels
epsCE  = [epsC(:),randn(size(epsC(:)))]*chol([1,rho_CE;rho_CE,1]); % account for correlation between consumption and systematic earnings shocks
epsE   = reshape(epsCE(:,2),size(epsC));

% simulate kappa panel
g_sim                 = mu_X(states) + sig_XS(states).*epsE + sig_Xid(states).*randn(Tsim,Ssim,Nsim); % [Tsim,Ssim,Nsim]
k_def_sim             = k_def(states); % [Tsim,Ssim]
k_iss_sim             = k_iss(states); % [Tsim,Ssim]
k_bar_sim             = k_bar(states); % [Tsim,Ssim]
i1                    = find(states(1,:)==1); % find periods corresponding to states 1 & 2
i2                    = find(states(1,:)==2);
k_sim                 = NaN(size(g_sim));
k_QML                 = NaN(size(g_sim));
last_issuance         = NaN(size(g_sim));
k_sim(1,i1,:)         = k_sim_init(i1_init(1:length(i1)),:);
k_sim(1,i2,:)         = k_sim_init(i2_init(1:length(i2)),:);
k_QML(1,i1,:)         = k_QML_init(i1_init(1:length(i1)),:);
k_QML(1,i2,:)         = k_QML_init(i2_init(1:length(i2)),:);
last_issuance(:,i1,:) = permute(repmat(last_issuance_init(i1_init(1:length(i1)),:),[1,1,Tsim]),[3,1,2]);
last_issuance(:,i2,:) = permute(repmat(last_issuance_init(i2_init(1:length(i2)),:),[1,1,Tsim]),[3,1,2]);
for t=1:Tsim-1
    
    k_sim(t+1,:,:)       = k_sim(t,:,:).*exp(g_sim(t,:,:));
    k_QML(t+1,:,:)       = k_QML(t,:,:);
    
    k_sim_temp           = k_sim(t+1,:,:); % temp: create 2D objects
    k_QML_temp           = k_QML(t+1,:,:);
    g_sim_temp           = g_sim(t,:,:);
    k_bar_temp           = k_bar(states(t,:))*ones(1,Nsim);
    states_temp          = states(t,:)'*ones(1,Nsim);
    
    default              = squeeze(k_sim(t,:,:)) <= k_def(states(t,:));
    k_sim_temp(default)  = k_bar_temp(default).*exp(g_sim_temp(default));
    k_QML_temp(default)  = k_bar_temp(default);
    issuance             = squeeze(k_sim(t,:,:)) >= k_iss(states(t,:));
    k_sim_temp(issuance) = k_bar_temp(issuance).*exp(g_sim_temp(issuance));
    k_QML_temp(issuance) = k_bar_temp(issuance);
    k_sim(t+1,:,:)       = k_sim_temp;
    k_QML(t+1,:,:)       = k_QML_temp;
    
    last_issuance_temp                = reshape(last_issuance(t:end,:,:),[Tsim-t+1,Ssim*Nsim]);
    last_issuance_temp(:,issuance(:)) = ones(Tsim-(t-1),1) * states_temp(issuance(:))';
    last_issuance(t:end,:,:)          = reshape(last_issuance_temp,[Tsim-t+1,Ssim,Nsim]);
    
end


default = k_sim <= repmat(k_def_sim,[1,1,Nsim]);


%--------------------------------------------------------------------------
% Table 5
%--------------------------------------------------------------------------


Fit     = griddedInterpolant({1:nS,k},lamDx);
DX_sim  = Fit(last_issuance,k_QML);
DX_sim(default) = NaN;

Fit     = griddedInterpolant({1:nS,k},lamSx);
SX_sim  = Fit(repmat(states,[1,1,Nsim]),k_sim);
SX_sim(default) = NaN;

Fit     = griddedInterpolant({1:nS,k},lamS);
S_sim   = Fit(repmat(states,[1,1,Nsim]),k_sim);
S_sim(default) = NaN;

QMLEV_sim = DX_sim./( DX_sim + SX_sim .* k_sim./k_QML);
QMLEV_sim(default) = NaN;

Fit     = griddedInterpolant({1:nS,k},lamDx);
DX_sim  = Fit(repmat(states,[1,1,Nsim]),k_sim); %%% true leverage
DX_sim(default) = NaN;

MLEV_sim = DX_sim./( DX_sim + SX_sim); %%% true leverage

Fit         = griddedInterpolant({1:nS,k},EP_R);
E_R_sim     = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k},VP_R');
V_R_sim    = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k},IV_ATM);
IV_sim      = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit        = griddedInterpolant({1:nS,k},IV_Skew);
Skew_sim   = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k,1:nS},CDS(:,:,:,60));
CDS60_sim   = Fit(repmat(states,[1,1,Nsim]),k_sim,last_issuance);

ER_sim      = E_R_sim - reshape(Rf(states,1),size(states));

Rcum_sim    = exp(log(S_sim(2:end,:,:))  - log(SX_sim(1:end-1,:,:)) + g_sim(1:end-1,:,:)) - 1;

R_MKT_sim   = trimmean(Rcum_sim,2,3); % expected returns are measured contemporaneously to mktcaps

data = [ QMLEV_sim(:) ER_sim(:) CDS60_sim(:) IV_sim(:) Skew_sim(:) reshape(repmat([NaN(1,Ssim); R_MKT_sim],[1,1,Nsim]),[Tsim*Ssim*Nsim,1]) V_R_sim(:)];
del = any(isnan(data),2);
data(del,:) = [];

mean_vec    = mean(data(:,1:5));
var_vec     = var(data(:,1:6));
var_vec(2)  = mean(data(:,7)) + var(data(:,2));
table_5     = [mean_vec, sqrt(var_vec)]';
G_m         = [mean_vec, var_vec];

tmp = MLEV_sim(:);
tmp(del) = [];
E_true_lev = mean(tmp);
SD_true_lev = std(tmp);


%--------------------------------------------------------------------------
% Table 6
%--------------------------------------------------------------------------


% Average time-series variances
ts_var = [mean(mean(squeeze(nanvar(QMLEV_sim)))), mean(mean(squeeze(nanmean(V_R_sim)+nanvar(ER_sim)))), ...
    mean(mean(squeeze(nanvar(CDS60_sim)))), mean(mean(squeeze(nanvar(IV_sim)))), mean(mean(squeeze(nanvar(Skew_sim))))]';

cs_var = NaN(5,1);

temp = QMLEV_sim;
temp(del) = NaN;
temp = nanvar(temp,[],3);
cs_var(1) = nanmean(temp(:));

temp = nan(size(QMLEV_sim));
temp(2:end,:,:) = Rcum_sim;
temp(del) = NaN;
temp(temp<prctile(temp(:),0.025)) = prctile(temp(:),0.025);
temp(temp>prctile(temp(:),99.975)) = prctile(temp(:),99.975);
temp = nanvar(temp,[],3);
cs_var(2) = nanmean(temp(:));

temp = CDS60_sim;
temp(del) = NaN;
temp = nanvar(temp,[],3);
cs_var(3) = nanmean(temp(:));

temp = IV_sim;
temp(del) = NaN;
temp = nanvar(temp,[],3);
cs_var(4) = nanmean(temp(:));

temp = Skew_sim;
temp(del) = NaN;
temp = nanvar(temp,[],3);
cs_var(5) = nanmean(temp(:));


%--------------------------------------------------------------------------
% Table 7
%--------------------------------------------------------------------------


Fit          = griddedInterpolant({1:nS,k},PD_q(:,:,60));
PD_Q_60_sim  = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit          = griddedInterpolant({1:nS,k},PD_p(:,:,60));
PD_P_60_sim  = Fit(repmat(states,[1,1,Nsim]),k_sim);

i1 = find(k>k_def(1),1);
i2 = find(k>k_def(2),1);
i3 = find(k>k_def(3),1);
k_small = [k(i1) k(i2) k(i3)]';
states_mat = repmat(states,[1,Nsim]);

Fit  = griddedInterpolant({1:nS,k,1:nS},LGD_q(:,:,:,1));
LGD_Q_1_sim = Fit(states_mat,k_small(states_mat),states_mat);

Fit  = griddedInterpolant({1:nS,k,1:nS},LGD_p(:,:,:,1));
LGD_P_1_sim = Fit(states_mat,k_small(states_mat),states_mat);

data = [PD_P_60_sim(:) PD_Q_60_sim(:) LGD_P_1_sim(:) LGD_Q_1_sim(:)];
data(del,:) = [];

mean_vec    = mean(data);
var_vec     = var(data(:,1:2));
table_7     = [mean_vec(1:2), sqrt(var_vec), mean_vec(3:4)]';


%--------------------------------------------------------------------------
% Table 8
%--------------------------------------------------------------------------


R_sim = ones(size(S_sim));
R_sim(2:end,:,:) = S_sim(2:end,:,:)./SX_sim(1:end-1,:,:);
Rsum  = zeros(size(S_sim));
sig_E = zeros(size(S_sim));
for t=12:Tsim
    Rsum(t,:,:) = exp(sum(log(R_sim(t-11:t,:,:))));
    sig_E(t,:,:) = std(log(R_sim(t-11:t,:,:)))*sqrt(12);
end
Rsum   = Rsum-1;
sig_D  = 0.05 + 0.25 * sig_E;
sig_V  = (1-QMLEV_sim) .* sig_E + QMLEV_sim .* sig_D;
DD_sim = (log(1./QMLEV_sim) + (Rsum - 1/2*sig_V.^2)*5)./(sig_V*sqrt(5));

Psim    = 5000;
IV_pct    = NaN(Tsim,Psim,10);
Skew_pct  = NaN(Tsim,Psim,10);
CDS_pct   = NaN(Tsim,Psim,10);

for t=1:Tsim
    
    IV_pct_i    = NaN(Psim,10);
    Skew_pct_i  = NaN(Psim,10);
    CDS_pct_i   = NaN(Psim,10);
    
    for jj=1:Psim
        X = zeros(1,10);
        
        Z = DD_sim(t,jj,:);
        bp = prctile(Z,10:10:90);
        
        Y = IV_sim(t,jj,:);
        X(1) = nanmean( Y( Z<=bp(1) ));
        for i=1:8
            X(i+1) = nanmean( Y( (Z>bp(i)) & (Z<=bp(i+1)) ));
        end
        X(10) = nanmean( Y( Z>bp(9) ));
        IV_pct_i(jj,:) = X;
        
        Y = Skew_sim(t,jj,:);
        X(1) = nanmean( Y( Z<=bp(1) ));
        for i=1:8
            X(i+1) = nanmean( Y( (Z>bp(i)) & (Z<=bp(i+1)) ));
        end
        X(10) = nanmean( Y( Z>bp(9) ));
        Skew_pct_i(jj,:) = X;
        
        Y = CDS60_sim(t,jj,:);
        X(1) = nanmean( Y( Z<=bp(1) ));
        for i=1:8
            X(i+1) = nanmean( Y( (Z>bp(i)) & (Z<=bp(i+1)) ));
        end
        X(10) = nanmean( Y( Z>bp(9) ));
        CDS_pct_i(jj,:) = X;
        
    end
    
    IV_pct(t,:,:) = IV_pct_i;
    Skew_pct(t,:,:) = Skew_pct_i;
    CDS_pct(t,:,:) = CDS_pct_i;
    
end

IV_pct = reshape(IV_pct,[Tsim*Psim,10]);
Skew_pct = reshape(Skew_pct,[Tsim*Psim,10]);
CDS_pct = reshape(CDS_pct,[Tsim*Psim,10]);

table_8 = [mean(IV_pct)', mean(Skew_pct)', 100*mean(CDS_pct)'];


